############### ICA ################
##### General component analysis ###
##### Author: Piatinskii M, 2019 ###
####################################

config.rodionov.L <- 7 # rodionov regime-shift minimum length (7 = regime shift can happens not ofter then 1 time in 7 years)
config.rodionov.alpha <- 0.05 # confidence level, alpha = 0.05 equals p = 0.95

source("rodionov.R")
# visualisation lib
library("factoextra")
library("gridExtra")
library("ggplot2")
library("reshape2")
library("rmarkdown")
library("dplyr")
library("foreach")
library("corrplot")
library("explor")
library("rgl")
library("NbClust")

# helper functions
plot.anomaly <- function(x, y, title) {
  r <- barplot(y, col = "black", space = 0.7, 
                  axis.lty = "solid", 
                  names.arg = x, axisnames = TRUE, 
                  ylab = "Anomalies", 
                  main = title)
  box(lty = "solid", col = "black")
  # creates a horizontal line at y = 0
  abline(h = 0)
}

# read input data and replace "," to "." and numeric type
data <- read.csv("input/data.csv")
# make ln(x + 1) linearization transformation
# for vegan or other packages that cannot scale=TRUE
data.log <- log(data+1)
data.log$year <- data$year

# step 1. Look up for anomaly diagnostics
data.anomaly <- data 
for (i in 2:ncol(data)) {
  mean <- mean(data.anomaly[,i])
  data.anomaly[,i] <- data.anomaly[,i] - mean(data.anomaly[,i])
}


png("output/anomaly/sprat.png", width=600, height=800, res=100, type="cairo")
par(mfrow = c(2,2))
plot.anomaly(data.anomaly$year, data.anomaly$sprat_ssb, "Sprat SSB")
plot.anomaly(data.anomaly$year, data.anomaly$sprat_rec, "Sprat Recruitment")
plot.anomaly(data.anomaly$year, data.anomaly$sprat_f, "Sprat Fishing mortality")
dev.off()

png("output/anomaly/env.png", width=600, height=800, res=100, type="cairo")
par(mfrow = c(2,2))
plot.anomaly(data.anomaly$year, data.anomaly$env_sst, "Sea surface temperature")
plot.anomaly(data.anomaly$year, data.anomaly$env_o2, "Disolved oxygen")
plot.anomaly(data.anomaly$year, data.anomaly$env_no3, "Nitrates (No3)")
plot.anomaly(data.anomaly$year, data.anomaly$env_po4, "Phosphates (Po4)")
dev.off()

png("output/anomaly/producers_consumers.png", width=800, height=800, res=100, type="cairo")
par(mfrow = c(2,3))
plot.anomaly(data.anomaly$year, data.anomaly$phyto_b, "Phytoplankton Biomass")
plot.anomaly(data.anomaly$year, data.anomaly$zoo_cold, "Zooplankton coldwater")
plot.anomaly(data.anomaly$year, data.anomaly$zoo_warm, "Zooplankton warmwater")
plot.anomaly(data.anomaly$year, data.anomaly$m_leidyi, "Mnemiopsis leidiy")
plot.anomaly(data.anomaly$year, data.anomaly$b_ovata, "Beroe ovata")
dev.off()

png("output/anomaly/predators.png", width=600, height=550, res=100, type="cairo")
par(mfrow = c(1,2))
plot.anomaly(data.anomaly$year, data.anomaly$turbot_crimea_b, "Turbot B")
plot.anomaly(data.anomaly$year, data.anomaly$bonito_b, "Bonito B")
dev.off()


# Performing PCA analysis using prcomp() default method
# Calculation done by singular value comomposition of the data matrix, not by using eigen on the covariance matrix.
# you can use "vegan" package and "rda" method to done eigen-value based approach
data.tmp <- data[,-1]
rownames(data.tmp) <- data$year

# fit PCA model
fit <- prcomp(data.tmp, center = TRUE, scale. = TRUE)

# get pca fit summary
fit.sum <- summary(fit)$importance

# get impact of each one in 1-3 PC
fit.impact <- as.data.frame(fit$rotation[,1:3])
fit.impact$data <- rownames(fit.impact)

# get weighted site score matrix for PC 1-3 in new coordinate system PCs
fit.score <- data.frame(year = as.numeric(data$year), pc1 = fit$x[,1], pc2 = fit$x[,2], pc3 = fit$x[,3])

# plot summary PCA images
# 1 - variance percentage by PCs
# Eigenvalues correspond to the amount of the variation explained by each principal component (PC)
# black-white pic
png("output/pca/bw/variance-parts.png", width=1900, height=1600, res=200, type="cairo")
print({
  fviz_eig(fit, barfill = "gray", barcolor = "black", choice = "variance")
})
dev.off()

# color pic
png("output/pca/col/variance-parts.png", width=1900, height=1600, res=200, type="cairo")
print({
  fviz_eig(fit, choice = "variance")
})
dev.off()

# 2 - variable graph. 
# Positive correlated variables point to the same side of the plot. Negative correlated variables point to opposite sides of the graph.
# black-white pic
png("output/pca/bw/variable-pc1-vs-pc2.png", width=1900, height=1600, res=200, type="cairo")
print({
  fviz_pca_var(fit) +  
    theme_minimal()
})
dev.off()

# color pic
png("output/pca/col/variable-pc1-vs-pc2.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_var(fit, col.var = "contrib", axes = c(1,2), gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel=TRUE)})
dev.off()

png("output/pca/col/variable-pc1-vs-pc3.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_var(fit, col.var = "contrib", axes = c(1,3),  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel=TRUE)})
dev.off()

png("output/pca/col/variable-pc2-vs-pc3.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_var(fit, col.var = "contrib", axes = c(2,3),  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel=TRUE)})
dev.off()

# 3 - contribution of variables by PC 1-3
# black-white pic
png("output/pca/bw/contrib-pc13.png", width=1900, height=1600, res=200, type="cairo")
print({grid.arrange(fviz_contrib(fit, choice = "var", axes = 1, color = "black", fill = "gray"), 
  fviz_contrib(fit, choice = "var", axes = 2, color = "black", fill = "gray"), 
  fviz_contrib(fit, choice = "var", axes = 3, color = "black", fill = "gray"),
  nrow=3)})
dev.off()

# colored pic
png("output/pca/col/contrib-pc13.png", width=1900, height=1600, res=200, type="cairo")
print({grid.arrange(fviz_contrib(fit, choice = "var", axes = 1), 
                    fviz_contrib(fit, choice = "var", axes = 2), 
                    fviz_contrib(fit, choice = "var", axes = 3),
                    nrow=3)})
dev.off()

# todo

png("output/pca/col/factor-contrib-year.png", width=1900, height=1600, res=200, type="cairo")
print({
  grid.arrange(fviz_contrib(fit, choice = "ind", axes = 1), 
               fviz_contrib(fit, choice = "ind", axes = 2), 
               fviz_contrib(fit, choice = "ind", axes = 3),
               nrow=3)
})
dev.off()

png("output/pca/bw/factor-contrib-year.png", width=1900, height=1600, res=200, type="cairo")
print({
  grid.arrange(fviz_contrib(fit, choice = "ind", axes = 1, color = "black", fill = "gray"), 
               fviz_contrib(fit, choice = "ind", axes = 2, color = "black", fill = "gray"), 
               fviz_contrib(fit, choice = "ind", axes = 3, color = "black", fill = "gray"),
               nrow=3)
})
dev.off()


# fviz_contrib(fit, choice = "var", axes = 1:3, color = "black", fill = "gray")

# 4 - variance and year historical changes
# black-white pic
png("output/pca/bw/biplot-year-variance.png", width=1900, height=1600, res=200, type="cairo")
#print({fviz_pca_biplot(fit, repel=TRUE, col.ind = "#696969")})
print({fviz_pca_biplot(fit, repel=TRUE, palette = "gray", col.ind = "gray", col.var = "black")})
dev.off()

# colored pic
png("output/pca/col/biplot-year-variance-pc1-2.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_biplot(fit, axes = c(1,2), repel=TRUE, col.ind = "#696969")})
dev.off()

png("output/pca/col/biplot-year-variance-pc1-3.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_biplot(fit, axes = c(1,3),repel=TRUE, col.ind = "#696969")})
dev.off()

png("output/pca/col/biplot-year-variance-pc2-3.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_pca_biplot(fit, axes = c(2,3),repel=TRUE, col.ind = "#696969")})
dev.off()

# 5 - final plot PC1-3 variance by year vector
m <- melt(fit.score, id="year")
# black white pic
png("output/pca/bw/pca-regime-shift.png", width=1900, height=1600, res=200, type="cairo")

print({ggplot(data=m, aes(x=year, y=value, colour=variable, size=variable, alpha=variable)) +
  geom_line(aes(linetype=variable)) + 
  geom_point(size = 2) + 
  scale_linetype_manual(values=c("solid", "dashed","dotted"),labels=c("PC1","PC2","PC3")) + 
  scale_color_manual(values=c("#000000", "#2e2e2e", "#4d4d4d"))+
  scale_size_manual(values=c(1.2, 1.0, 0.7),labels=c("PC1", "PC2", "PC3")) + 
  scale_alpha_manual(values=c(0.9,0.8,0.7)) + 
  xlab("Year") + 
  ylab("PCA variance") +
  #annotate("rect", xmin=2007, xmax=2011, ymin=-Inf, ymax=Inf, alpha=0.1, fill="blue") + 
  geom_hline(yintercept = 0)})

dev.off()

# color pic
png("output/pca/col/pca-regime-shift.png", width=1900, height=1600, res=200, type="cairo")

print({ggplot(data=m, aes(x=year, y=value, colour=variable, size=variable, alpha=variable)) +
    geom_line(aes(linetype=variable)) + 
    geom_point(size = 2) + 
    scale_color_manual(values=c("red", "blue", "darkgreen"))+
    scale_size_manual(values=c(1.1, 1.0, 0.9),labels=c("PC1", "PC2", "PC3")) + 
    scale_alpha_manual(values=c(1,1,1)) + 
    xlab("Year") + 
    ylab("PCA variance") +
    #annotate("rect", xmin=2007, xmax=2011, ymin=-Inf, ymax=Inf, alpha=0.1, fill="blue") + 
    geom_hline(yintercept = 0)})

dev.off()

# 6 - process "regime shift detection" by rodionov, 2004 algo
regime.shift <- rstars(data.timeseries = fit.score, pValue = config.rodionov.alpha, l.cutoff = config.rodionov.L)
rod <- list(pc1 = c(), pc2 = c(), pc3 = c(), year = fit.score$year)

# rod$pc1 <- rodionov(fit.score$PC1, L = config.rodionov.L , p = config.rodionov.alpha)
# rod$pc2 <- rodionov(fit.score$PC2, L = config.rodionov.L , p = config.rodionov.alpha)
# rod$pc3 <- rodionov(fit.score$PC3, L = config.rodionov.L , p = config.rodionov.alpha) # fail performing at L=10

rod.df <- as.data.frame(regime.shift$mean)
#names(rod.df) <- c("pc1-rod", "pc2-rod", "pc3-rod", "year")

m <- melt(rod.df, id="year")
# black white pic
png("output/pca/bw/pca-regime-shift-radeonov.png", width=1900, height=1600, res=200, type="cairo")

print({ggplot(data=m, aes(x = year, y=value, size=variable, alpha=variable, linetype=variable)) + 
  geom_line() + 
  scale_linetype_manual(values=c("solid", "dashed","dotted")) + 
  scale_color_manual(values=c("#000000", "#2e2e2e", "#4d4d4d"))+
  scale_size_manual(values=c(1.2, 1.0, 0.7)) + 
  scale_alpha_manual(values=c(0.9,0.8,0.7)) + 
  geom_hline(yintercept = 0) +
  geom_point() + 
  xlab("Year") + 
  ylab("Rodionov PCs values")})

dev.off()

# color pic
png("output/pca/col/pca-regime-shift-radeonov.png", width=1900, height=1600, res=200, type="cairo")

print({ggplot(data=m, aes(x = year, y=value, size=variable, alpha=variable, colour = variable)) + 
    geom_line() + 
    scale_color_manual(values=c("red", "blue", "darkgreen"))+
    scale_size_manual(values=c(1.1, 1.0, 0.9)) + 
    scale_alpha_manual(values=c(1,1,1)) + 
    geom_hline(yintercept = 0) +
    geom_point() + 
    xlab("Year") + 
    ylab("Rodionov PCs values")})

dev.off()

# 7 - final plot with regime shift
merged <- left_join(fit.score, rod.df, by = "year")
m <- melt(merged, id="year")
png("output/pca/col/pca-regime-shift-data-radeonov.png", width=1900, height=1600, res=200, type="cairo")

print({ggplot(data=m, aes(x=year, y=value, color=variable, size=variable, alpha=variable)) +
  geom_line() + 
  scale_linetype_manual(values=c("solid", "solid", "solid", "dashed","dashed", "dashed")) +
  scale_colour_manual(values=c("#e84a3f", "#3fbee8", "#479156", "#e84a3f", "#3fbee8", "#479156")) +
  scale_size_manual(values=c(1.4, 1.0, 0.8, 0.4, 0.4, 0.4)) + 
  scale_alpha_manual(values=c(1,1,1, 0.7,0.7, 0.7)) + 
  xlab("Year") + 
  ylab("PCA variance") +
  #annotate("rect", xmin=2007, xmax=2011, ymin=-Inf, ymax=Inf, alpha=0.1, fill="blue") + 
  geom_hline(yintercept = 0)})

dev.off()

### Perform K-means clustering procedure
# find optimum number of clusters
png("output/pca/col/pca-kmeans-cluster-silhouette.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_nbclust(fit$x, kmeans, method = "silhouette")})
dev.off()

# compute cluster optimum by NbClust()
nb <- NbClust(fit$x, distance = "euclidean", min.nc = 2, max.nc = 5, method = "kmeans")
# set numbers there 
clust.numb <- 3
# or use
# clust.numb <- length(unique(nb$Best.partition))

kc12 <- kmeans(fit$x[,c(1,2)], clust.numb)
kc13 <- kmeans(fit$x[,c(1,3)], clust.numb)
kc23 <- kmeans(fit$x[,c(2,3)], clust.numb)

# TODO: The number of clusters can be found automaticly by eclust() method
ec12 <- eclust(fit$x[, c(1,2)], "kmeans", nboot = 1000)
ec13 <- eclust(fit$x[, c(1,3)], "kmeans", nboot = 1000)
ec23 <- eclust(fit$x[, c(2,3)], "kmeans", nboot = 1000)

# plot(fit$x[,c(1,2)], col=factor(kc12$clust), main = "Kmeans clustering", pch=16)
# abline(h = 0, lty = "dashed")
# abline(v = 0, lty = "dashed")
# text(fit$x[,1], fit$x[,2]+0.2, rownames(fit$x), cex = 0.9)
# the same in fviz
png("output/pca/col/pca-kmeans-pc-1-2.png", width=1900, height=1600, res=200, type="cairo")
print({fviz_cluster(kc12, data = fit$x[,c(1,2)], ggtheme = theme_minimal()) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed")})
dev.off()

png("output/pca/col/pca-kmeans-pc-1-3.png", width=1900, height=1600, res=200, type="cairo")
fviz_cluster(kc13, data = fit$x[,c(1,3)], ggtheme = theme_minimal()) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed")
dev.off()

png("output/pca/col/pca-kmeans-pc-2-3.png", width=1900, height=1600, res=200, type="cairo")
fviz_cluster(kc23, data = fit$x[,c(2,3)], ggtheme = theme_minimal()) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed")
dev.off()

# 8 - try to build cluster dendrogram with same years by PC variations
hc <- hclust(dist(cbind(fit$x[,1], fit$x[,2], fit$x[,3])), method = "ward.D2")
png("output/pca/bw/pca-cluster-dendrogram.png", width=1900, height = 1600, res = 200, type="cairo")
plot(hc, main = "PCA cluster dendrogram by PC variations impact", xlab = "Obs. year", ylab = "Cluster height")
rect.hclust(hc, k=clust.numb, border="gray")
dev.off()

# all input data correlation test and p-value signif.level processing
png("output/correlation-test-bw.png", width=1900, height = 1600, res = 200, type="cairo")
cor.mat <- cor(data[,-1])
p.mat <- cor.mtest(data[,-1])$p
col <- colorRampPalette(gray.colors(3, start = 0.5, end = 0.9))
corrplot(cor.mat, method = "color", col = gray.colors(200),
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90,
         p.mat = p.mat, sig.level = 0.05, insig = "blank", 
         diag = FALSE)
dev.off()

png("output/correlation-test-col.png", width=1900, height = 1600, res = 200, type="cairo")
cor.mat <- cor(data[,-1])
p.mat <- cor.mtest(data[,-1])$p
corrplot(cor.mat, method = "color",
         type = "upper", order = "hclust", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90,
         p.mat = p.mat, sig.level = 0.05, insig = "blank", 
         diag = FALSE)
dev.off()

# 9 - traffic-light plot 

# get data dimensions
d <- dim(data)
data.names <- names(data)[-1]
# get dataset without year column
data.noyear <- data[,-1]

# prepare quantile statistics matrix
quant <- matrix(nrow = 5, ncol = (ncol(data.noyear)))
# fill matrix with quantile values for levels: 0.2, 0.4, 0.6, 0.8, 1
for (i in 1:ncol(data.noyear)) {
  quant[,i] <- quantile(data.noyear[,i], probs = c(0.2, 0.4, 0.6, 0.8, 1), na.rm = T)
}

# create quantile-relation matrix categorization
# focal point: create factor matrix with determination 
# to quantile intervals 
# ex, 1 means that value is in 0 ... 0.2 quantile of full factor row
quant.matrix <- data.noyear
for (i in 1:(nrow(quant.matrix))) {
  for (j in 1:(ncol(quant.matrix))) {
    cel <- quant.matrix[i, j]
    if (is.na(cel))
      quant.matrix[i,j] <- TRUE
    else if (cel < quant[1, j]) # cel in 1st quantile value between 0 ... 0.2 in col
      quant.matrix[i, j] <- 1 # set 1st quantile to output
    else if (cel >= quant[1, j] && cel < quant[2, j]) 
      quant.matrix[i, j] <- 2
    else if (cel >= quant[2, j] && cel < quant[3, j]) 
      quant.matrix[i, j] <- 3
    else if (cel >= quant[3, j] && cel < quant[4, j]) 
      quant.matrix[i, j] <- 4
    else if (cel >= quant[4, j])
      quant.matrix[i, j] <- 5
    
  }
}

# get 1st PC values
pc1 <- fit.impact[,1]
# create ascending index from PC1 values (from negative to positive)
pc1.idx <- order(pc1)

# order quantile-transformed data columns by pc1 impact index
quant.matrix.order <- quant.matrix[pc1.idx]
# create matrix object from ordered quantiles by pc1 impact
matrix <- as.matrix(quant.matrix.order)

png("output/pca/col/traffic-light.png", width=1900, height=1600, res=200, type="cairo")
# finally prepare pics
x <- 1:(nrow(data.noyear))
y <- 1:(ncol(data.noyear))
op <- par(mar = c(3, 3, 3, 6), oma = c(0.5,2,0,0), xpd = TRUE)
image(x, y, z = matrix, zlim = c(1,5), 
      col = c("green", "yellowgreen", "yellow", "gold", "red"), 
      axes = FALSE, xlab = "Year", ylab = "")
axis(1, at = seq(1, nrow(data.noyear), by = 1), tick = FALSE, labels = data$year, 
     cex.axis = 0.7, las = 3)
axis(2, at = seq(1, ncol(data.noyear), by = 1), tick = FALSE, 
     labels = data.names[pc1.idx],
     cex.axis = 0.8, las = 1, padj = 1) 
box()
title(main = "Traffic Light Plot", font.main = 2)
legend(x = max(x)+0.5, y = max(y), title = "Quantile colors:",
       legend = c("0 ... 0.2", "0.2 ... 0.4", "0.4 ... 0.6", "0.6 ... 0.8", "0.8 ... 1"),
       fill = c("green", "yellowgreen", "yellow", "gold", "red"),
       border =  c("green", "yellowgreen", "yellow", "gold", "red"), 
       cex = 0.7, bty = "n")
par(op)
dev.off()



# render output reports for black-white and color scheme
#render("Report_bw.Rmd")
render("Report_col.Rmd")

#corrplot(cor.results, method="number")
#explor(fit)